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Elastic effects in a model of disordered nematic elastomers are numerically investigated in two 
dimensions. Networks crosslinked in the isotropic phase exhibit unusual soft mechanical response 
against stretching. It arises from gradual alignment of orient at ionally correlated regions that are 
elongated along the director. A sharp crossover to a macroscopically aligned state is obtained on 
further stretching. The effect of random internal stress is also discussed. 

PACS numbers: 61.30.Cz, 61.41. +e, 64.70.Md 
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Nematic elastomers and gels exhibit rich mechanical 
effects due to elasticity-orientation coupling [QJ2). While 
a considerable number of theoretical studies has been di- 
rected to homogeneous systems, nematic elastomers are 
often in a highly non-uniform polydomain state, in which 
the correlation length for the director orientation is typ- 
ically of micron scales. Polydomain networks show un- 
usual non-linear clastic response against stretching, of- 
ten with an extremely low stress over a large interval of 
strain ||^,||,§,0J§1. As the strain is increased, the direc- 
tors gradually rotate toward the direction of stretching 
until a macroscopically aligned state is attained. This 
structural change is called the polydomain-monodomain 
(P-M) transition. Attempting to describe the presum- 
ably equilibrium polydomain textures, Terentjev and 
coworkers ||,[lC]|llJ proposed a random-field model anal- 
ogous to those for random anisotropy magnets. They ar- 
gued that crosslinkers of anisotropic shapes act as sources 
of quenched disorder. On the other hand, the mechan- 
ical response is not yet well understood. It is known 
that elasticity-mediated long-range interactions among 
spatial inhomogeneities are crucial in systems such as 
metallic alloys [|l 2[Jl 3f| and gels juj . For polydomain net- 
works, the role of elastic interactions among orientation- 
ally correlated regions ("domains") is yet to be clarified. 
In this Rapid Communication, we numerically investi- 
gate the mechanical response and the domain structure of 
model nematic networks incorporating both rubber elas- 
ticity and quenched random anisotropy. Unusual soft re- 
sponse is obtained and is explained in terms of the elastic 
interaction. We briefly discuss the effect of random in- 
ternal stress as another kind of quenched disorder which 
can destroy long-range orientational order Q . 

The total free energy of our model system is of the 
form F = F e i + Fr + Fp, where F e i, Fr, and Fp are, 
respectively, the rubber-elastic, random disorder, and 
Frank contributions. We assume networks brought deep 
into the nematic phase after crosslinking in the isotropic 
phase, and apply the affine-deformation theory of ne- 
matic rubber elasticity due to Warner et.al. |£6[. Then 
F e i is written in terms of the symmetric deformation ten- 
sor Wij = {dr i /dr k v ){drj/dr k ) ), where r° and r-j are the 
Cartesian coordinates of the material point at the mo- 
ment of crosslinking and after deformation, respectively. 



Summation over repeated indices is implied throughout 
the paper unless otherwise stated. It is convenient to 
rewrite the original form fl6[| of F e i using the tensor 
Qij — TiiUj — Sij/d, where n is the director, to obtain M] 



*rf = f / driWu-aQijWij). 



(1) 



The dimensionless coupling constant a(> 0) is deter- 
mined by chain anisotropy and does not exceed d/(d—l). 
The modulus /i is given by fcgT multiplied by the 
crosslink number density and a numerical prefactor (~ 1) 
which is weakly dependent on the temperature. We con- 
sider the incompressible limit and impose the constraint 
detW = 1. The disorder free energy is assumed in the 
form given in |^,|l^,^l) , and is rewritten as 



Fr = / dr P l3 Q 



(2) 



where Pij is a symmetric, traceless, Gaussian random 
tensor with vanishing quenched average ({Pij(r)) = 0) 
and with variance 



(3) 



(P ij (q)P kl (-q)) = ufs ik Sji +5u5 jk - 
For the Frank free energy we assume the form 

Ff = f J dr (Vn)2 - (4) 

Here we treat the two-dimensional case for numeri- 
cal and analytical advantages. Then, in the absence of 
elasticity, our model reduces to the random-anisotropy 
2D XY model by regarding the unit vector m = 
(2Q xx ,2Q xy ) = (cos 26, sin 29) as the spin variable, where 
9 is the director orientation defined by n = (cos 9, sin ff). 
We consider deformations of the form = Air? + m (no 
summation) where Aa; = A and X y — 1/A express the 
average deformation, and u = u(r) denotes the internal 
displacement. Cooling into the nematic phase tends to 
induce spontaneous elongation along the director Jl],|| . If 
the directors are uniformly aligned along the £-axis, the 
elastic free energy ([!]) is minimized at A = A m and u = 
with 
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A, 



l + a/2 
1 - a/2 



1/4 
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However, if the random field is sufficiently strong, the 
system under no external stress energetically favors a 
macroscopically isotropic state with A = 1. This can 
be identified with the polydomain state. Our questions 
concern how domains spontaneously deform and to what 
degree the elastic free energy is reduced in this highly 
non-uniform state. 

The mechanical response was numerically simulated by 
varying the macroscopic strain A and minimizing the free 
energy with respect to n and u for each value of A. We 
solved the Langevin equation for the director, 



dn 
~dt 



(J — nn) 
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where n is an uncorrelated Gaussian thermal noise in- 
troduced to facilitate structural evolution. Without the 
noise the minimization process would stop at one of the 
local minima close to the initial configuration. After ap- 
proaching the global minimum we turned off the noise as 
explained below. The displacement u was determined by 
solving the non-linear equation S(F e i + F v )/Su = with 
a relaxation method, where F v is an artificial free energy 
functional of it, which penalizes volume change. In this 
way the local volume was kept constant with errors below 
1% throughout the runs. Periodic boundary conditions 
were imposed on n and u. The simulation was performed 
on a 128 2 square lattice with the grid size Ax = 1. We 
set K = 4 and U = 1 for all the runs, whereas \i and a 
were varied for different runs. In each run the external 
strain A was slowly increased after an initial equilibration 
stage at A = 1. Occasionally, we stopped the increase 
of A and turned off the thermal noise for an interval of 
time. Thus a single run consisted of alternating periods 
of annealing (with increasing strain) and quenching. In 
each quench period we computed the spatially-averaged 
free energy density f = f e i + fn, + If and the orienta- 
tion S — (2Q XX ) — (003 26*). This procedure enabled us 
to approximately minimize the free energy at numerous 
values of A in reasonable computational time. For a fur- 
ther check, we then decreased A back from a large value 
in a similar manner. A small hysteresis was obtained but 
it does not affect the description below. 

In Fig.l we show the strain-stress and strain- 
orientation relations for several values of a with /j,a 2 = 4 
fixed. In both plots we can see a sharp crossover around 
A = A m (a). Below A m the total stress df/dX is vanish- 
ingly small and slightly positive. The average orientation 
increases almost linearly in the same region. The free en- 
ergy densities are plotted in Fig. 2. The elastic free energy 
has a slightly negative slope in the region A < A m , while 
the disorder free energy has a positive slope and makes 
the total stress slightly positive. We chose the parame- 
ters so that the Frank contribution is much smaller than 
/ia 2 , which is considered to be the case in typical ex- 
periments. We also studied a few cases with stronger or 



weaker elastic effects. For larger values of \io? the shapes 
of the strain-elastic stress and strain-orientation curves 
were almost unchanged. For cases with \io? < 0.2 these 
two curves exhibited less sharp crossovers. 

In order to understand the origin of the soft response 
it is useful to examine the structure of the polydomain 
state at A = 1, for which an analytical treatment is 
possible in the weak coupling case a < 1. We expand 
AF e i — F e i [u] — F e i [0] with respect to Vit to obtain 



AF el =n dr 



1 / dui du 



4 V dr .; dr 



dui 



dri 



(7) 



Eliminating the elastic field using the conditions of me- 
chanical equilibrium SAF e i /Su = and incompressibility 
V • u = 0, we have a non-local elastic interaction among 
orientational inhomogeneities. We define new variables 
Qi{r) and Q2{t) through their Fourier transforms, 

Qi(q) = sbx(2(p)Q xx (q) - cos{2ip)Q xy {q) 1 (8) 
Q 2 (<?) = cos{2ip)Q xx (q) + sin(2ip)Q xy (q), (9) 

where <p is the azimuthal angle of the wave- vector q = 
q(cos<p,smip). Then the average free energy density 
reads lOl 



/dU=i =/*(!- y (Qi) 



(10) 



to order a 2 . Note that Qi and Q2 satisfy (Q\ + Q|) = 
(Qlx + Ql y ) = 1/4. We have (Qf) = (Q 2 ) = 1/8 in the 
absence of the elastic coupling. In its presence, asym- 
metry (Q 2 ) > (Q 2 ) arises to reduce the elastic free en- 
ergy (lOh. In the elasticity-dominated limit where /ia 2 
is much larger than the disorder and the Frank free en- 
ergy densities, we expect (Q 2 ) — > 1/4, (Q|) ~> 0, and 
fei\x=i — > — a2 /8). Indeed these are numerically 
confirmed as shown in the next paragraph. On the other 
hand, the elastic free energy density under the uniform 
deformation with A = X m is also given by /i(l — a 2 /8) to 
order a 2 . Thus, in the above limit, the P-M transition 
accompanies only a small change of order a 3 in the elas- 
tic free energy. To see how each domain is deformed at 
A = 1, we consider the local elastic stress which is given 



as <7; 



djUi 



xQij) from (0). After some 



calculation, its variance in the mechanical equilibrium is 
obtained as 



<^- 2 ) = 2^0? {Ql) 



(11) 



In the elasticity-dominated limit the variance of the 



quantity /i Oj 



diUj 



aQij vanishes due to 



the factor (Q 2 ) in ([il|), which means that each part of 
the system is elongated by 1 + a/A (~ A m ) times along 
the local director. This, together with the numerical re- 
sult on the mechanical response, supports the following 
simple picture : In the polydomain state each domain is 
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uniaxially elongated by A m times along the local director, 
and thus the elastic free energy is equal to that for the 
monodomain state at A = A m (Fig. 3). The P-M transi- 
tion in the region 1 < A < A m proceeds via rotation of 
domains and does not change the elastic free energy. 

Next we present numerical results on the polydomain 
structure at A = 1, which was studied through the cor- 
relation function G(r) = 2 (Qij(r)Qij(Q)) and the degree 
of structural asymmetry A = (Qf) — (<3|)- To acceler- 
ate computation of the elastic field we assumed a weak 
coupling a = 0.1 and solved 5AF e [/5u — under the 
constraint V • u = using FFT, instead of the relax- 
ation method above. The amplitude of the thermal noise 
was set constant in an initial stage and then gradually 
reduced to zero at a constant rate. The correlation func- 
tion is computed for the final state and averaged over 20 
independent runs for each set of parameters. Runs were 
sufficiently long to insure that the initial configurations 
with uniform and random orientations give indistinguish- 
able results for G(r). Shown in Figs. 4 and 5 are the cor- 
relation function and the correlation length R defined by 
G(R)/G(0) — 1/2. The elastic coupling increases the cor- 
relation length without qualitatively affecting the form of 
the correlation function. We could not deduce a quantita- 
tive decay law for G(r) from the relatively small number 
of samples, but the decay was slightly faster than expo- 
nential near the origin. For the non-elastic case the same 
feature was obtained in the Monte Carlo simulation by 
Gingras and Huse in the presence of thermal noise, 
while Yu et.al. jllj obtained exponential decay using free 
boundary conditions. Another important factor affecting 
G(r) is the disorder strength. More systematic study of 
the decay law is left to future work. In Fig. 5 the degree of 
asymmetry A is also shown. With increasing the magni- 
tude of the elastic interaction it approaches to the upper 
limit 1/4 as expected. 

Finally we discuss the effect of random internal stress 
arising from microscopic heterogeneities in the network 
structure, which are intrinsic to gels Jl9| . We restrict our 
discussion to the case A = 1 with small internal defor- 
mations. In the expansion of the clastic free energy with 
respect to Vm there will arise an additional term, 

AF el>R = J dr R ij (12) 

where Rij is the Gaussian random stress with (Rij(r)) = 
and 

(R ij (q)Rki(-q)) = ViSySki + V 2 {8 lk 8 j i + 5 a 8 jk ). (13) 

Eliminating the elastic field from AF e i + AF e | jj we have 
a new interaction free energy a J drRfQi, where Rf is 



defined using the shear component Rfj = Rij — RkkSij/d 



by an equation parallel to (|8|) as 
Rf(q) = S in(2<p)R s xx (q) 



(14) 



cos(2<p) R b xy (q). 

Treating this interaction as a weak perturbation as 
in pC| ], we can see that it renders the equilibrium cor- 
relation length finite even in the absence of the disorder 



free energy (3). We mention that Golubovic and Luben- 
sky jl5| discussed another mechanism of long-range- 
orientational-order breaking due to random stress. Their 
argument is based on the observation that the amplitude 
of thermal fluctuations around a uniformly aligned state 
diverges. Its relevance to the present case of nematic 
networks is limited in that their free energy does not ex- 
plicitly include the orientational degree of freedom. 

To summarize, we have numerically obtained a soft 
mechanical reponse during the P-M transition. It orig- 
inates from structural self-organization of domains due 
to the long-range elastic interaction, and should be dis- 
tinguished from the soft elasticity |p|,pT| of uniformly ori- 
ented networks. The elastic contribution to the stress is 
slightly negative in the transition region. We have found 
a positive disorder contribution to the stress. The elas- 
tic interaction is found to increase the correlation length. 
We have demonstrated that random internal stress acts 
as a random field on the director. Further experimen- 
tal and theoretical studies are necessary to examine its 
relevance to real polydomain textures. 
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strain X 

FIG. 1. Top: dimensionless total stress fi~ 1 df/d\. Bottom: orientation S = (2Q XX ) = (cos 28). Cases with different coupling 
strengths a = 0.2, 0.4, 0.6, 0.8 from left to right are compared with fia 2 — 4 fixed. The arrows indicate the corresponding values 
of A m . 
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FIG. 2. Free energy densities for fia 2 — 4 and a — 0.4. The total free energy in the polydomain regime has a positive but 
small slope due to the disorder contribution. The value of n is subtracted from the elastic free energy density. 
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FIG. 3. Schematic illustration of the P-M transition. The ellipses represent domains under spontaneous deformations (from 
circles at the moment of crosslinking) and the arrows in them indicate the local director orientations. Transition from polydo- 
main at A = 1 (left) to monodomain at A = A m (right) does not change the elastic free energy if every domain is elongated by 
A m times along the local director. 
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FIG. 5. Top: correlation length R. Bottom: structural asymmetry A = (Qi) — (QV)- 



5 



